% fields H for rho > a & rho < a

function [Hrho, Hphi, Hz, Erho, Ephi, Ez] = func_HE_inside(r, m, CONSTS, coef, p1)

    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eta   = CONSTS.eta;
    eps   = CONSTS.eps;
    g     = CONSTS.g;
    
%     if (isnan(p1))
%         Hrho = 0;
%         Hphi = 0;
%         Hz   = 0;
%         Erho = 0;
%         Ephi = 0;
%         Ez   = 0;
%         return;
%     end

    
    B1 = coef(1);
    B2 = coef(2);

    [q1 q2] = func_qk_from_p(p1, CONSTS);
    [n1 n2] = func_nk_from_p(p1, CONSTS, false);
    [alfa_1 alfa_2 beta_1 beta_2] = func_alfa_beta_from_p(p1, CONSTS, false);

    Q1 = k0*a*q1;
    Q2 = k0*a*q2;
    bessel_m_q1 = besselmx(double('J'), m, Q1*r, 0); %besselj(m, Q1*r);
    bessel_mp1_q1 = besselmx(double('J'), m+1, Q1*r, 0); %besselj(m+1, Q1*r);
    bessel_m_q2 = besselmx(double('J'), m, Q2*r, 0); %besselj(m, Q2*r);
    bessel_mp1_q2 = besselmx(double('J'), m+1, Q2*r, 0); %besselj(m+1, Q2*r);
    
    Hrho  = -1i*(B1*(p1*bessel_mp1_q1 - n1*beta_1*m*(bessel_m_q1./(Q1*r))) + ...
                         B2*(p1*bessel_mp1_q2 - n2*beta_2*m*(bessel_m_q2./(Q2*r))));
    
    Hphi  = -(B1*n1*(bessel_mp1_q1 - beta_1*m*(bessel_m_q1./(Q1*r))) + ...
                      B2*n2*(bessel_mp1_q2 - beta_2*m*(bessel_m_q2./(Q2*r))));
    
    Hz    = -(B1*q1*bessel_m_q1 + B2*q2*bessel_m_q2);
    

    Erho  = -(B1*(((n1*p1+g)/eps)*bessel_mp1_q1 + alfa_1*m*(bessel_m_q1./(Q1*r))) + ...
                      B2*(((n2*p1+g)/eps)*bessel_mp1_q2 + alfa_2*m*(bessel_m_q2./(Q2*r))));
    
    Ephi  = 1i*(B1*(bessel_mp1_q1 + alfa_1*m*(bessel_m_q1./(Q1*r))) + ...
                        B2*(bessel_mp1_q2 + alfa_2*m*(bessel_m_q2./(Q2*r))));
    
    Ez    = (1i/eta)*(B1*n1*q1*bessel_m_q1 + B2*n2*q2*bessel_m_q2);
    
end
